#include "saveGmsh.h"

using namespace std;

//------------------------------------------------------------------------------
SaveOutputGmsh::SaveOutputGmsh(Mesh* mymesh,FsContainer* myfs){
  mesh = mymesh;
  fs   = myfs;
}
//------------------------------------------------------------------------------
void SaveOutputGmsh::saveOutput(const SolContainer* sol,int Step,double t,int nbVar,string Name){

  int nbComps = 1;
  switch (nbVar) {
    case 1: nbComps = 1;break;
    case 2: nbComps = 3;break;  
    default:break;
  }
  
  stringstream tmp; tmp << Step;

  string filename = "output/"+Name+".t"+tmp.str()+".msh";
  FILE* out = fopen(filename.c_str(),"w");

  // --- Write mesh format ---

  int one = 1;
  fprintf(out,"$MeshFormat \n2.2 1 8 \n");
	fwrite(&one, sizeof(int), 1,out);			
	fprintf(out,"$EndMeshFormat \n");
    
  // --- Write interpolation scheme ---
  
  for(FsContainer::const_iterator fsIt = fs->begin();fsIt!=fs->end();++fsIt){
    
    int nbENP          = (*fsIt)->getNbENP();
    const DGMatrix* pp = (*fsIt)->getEpp();
    
    fprintf(out,"$InterpolationScheme\n");
    fprintf(out,"\"DGInterpolationScheme\"\n");
    fprintf(out,"1\n");
    fprintf(out,"%i\n",3);                                                      // gmsh elmt type   //NOT YET IMPLEMENTED for quad
    fprintf(out,"2\n");                                                         // nb of matrices
    fprintf(out,"%i %i\n",nbENP,nbENP);                                         // size of interp DGMatrix
    for(int i=0;i<nbENP;i++){
      for(int j=0;j<nbENP;j++) fprintf(out,"%g ",(*fsIt)->getERefNodesCoef(i,j)); 
      fprintf(out,"\n"); 
    }
    fprintf(out,"%i 2\n",nbENP);                                                // fs size & dim
    for(int i=0;i<nbENP;i++) fprintf(out,"%i %i\n",(int)(*pp)(i,0),(int)(*pp)(i,1)); 
    fprintf(out,"\"$EndInterpolationScheme\"\n");
    
  }
  
  // --- Write solution ---
        
  int eT=0;
  for(FsContainer::const_iterator fsIt = fs->begin();fsIt!=fs->end();++fsIt,eT++){
    
    // --- Write infos data ---
    
    fprintf(out,"$ElementNodeData\n");
    fprintf(out,"2\n");
    fprintf(out,"\"%s \"\n",Name.c_str());
    fprintf(out,"\"DGInterpolationScheme\"\n");
    fprintf(out,"1\n");
    fprintf(out,"%g\n",t);
    fprintf(out,"4\n");
    fprintf(out,"%i\n",Step);
    fprintf(out,"%i\n",nbComps);                                                // nbEq //NOT YET IMPLEMENTED nb equations != 1
    fprintf(out,"%i\n",mesh->getNbElmt(eT));
    fprintf(out,"0\n");
    
    // --- Write solution ---
    
    int nbENP = (*fsIt)->getNbENP();
 
    for(int eID=0;eID<mesh->getNbElmt(eT);eID++){
      int elmtGmshID = mesh->getGmshTag(eT,eID);
      fwrite(&elmtGmshID,sizeof(int),1,out);
      fwrite(&nbENP,sizeof(int),1,out);
      for(int k=0;k<nbENP;k++){
        for(int dim=0;dim<nbVar;dim++){
          double u = (*sol)[eT](k,eID*nbVar+dim);
          fwrite(&u,sizeof(double),1,out);
        }
        if(nbComps==3){double u = 0.; fwrite(&u,sizeof(double),1,out);}           // add 0 because vector is 3D!
      }
    }

    fprintf(out,"\"$EndElementNodeData\"\n");

  }  
  
  fclose(out);
}
//------------------------------------------------------------------------------
void SaveOutputGmsh::saveOutput(const DGMatrix* sol,int Step,double t,int nbVar,string Name){
  
  int nbComps = 1;
  switch (nbVar) {
    case 1: nbComps = 1;break;
    case 2: nbComps = 3;break;  
    default:break;
  }
  
  stringstream tmp; tmp << Step;
  
  string filename = "output/"+Name+".t"+tmp.str()+".msh";
  FILE* out = fopen(filename.c_str(),"w");
  
  // --- Write mesh format ---
  
  int one = 1;
  fprintf(out,"$MeshFormat \n2.2 1 8 \n");
	fwrite(&one, sizeof(int), 1,out);			
	fprintf(out,"$EndMeshFormat \n");
  
  // --- Write interpolation scheme ---
  
  for(FsContainer::const_iterator fsIt = fs->begin();fsIt!=fs->end();++fsIt){
    
    int nbENP          = (*fsIt)->getNbENP();
    const DGMatrix* pp = (*fsIt)->getEpp();
    
    fprintf(out,"$InterpolationScheme\n");
    fprintf(out,"\"DGInterpolationScheme\"\n");
    fprintf(out,"1\n");
    fprintf(out,"%i\n",3);                                                      // gmsh elmt type   //NOT YET IMPLEMENTED for quad
    fprintf(out,"2\n");                                                         // nb of matrices
    fprintf(out,"%i %i\n",nbENP,nbENP);                                         // size of interp DGMatrix
    for(int i=0;i<nbENP;i++){
      for(int j=0;j<nbENP;j++) fprintf(out,"%g ",(*fsIt)->getERefNodesCoef(i,j)); 
      fprintf(out,"\n"); 
    }
    fprintf(out,"%i 2\n",nbENP);                                                // fs size & dim
    for(int i=0;i<nbENP;i++) fprintf(out,"%i %i\n",(int)(*pp)(i,0),(int)(*pp)(i,1)); 
    fprintf(out,"\"$EndInterpolationScheme\"\n");
    
  }
  
  // --- Write solution ---
  
  int eT=0;
  for(FsContainer::const_iterator fsIt = fs->begin();fsIt!=fs->end();++fsIt,eT++){
    
    // --- Write infos data ---
    
    fprintf(out,"$ElementNodeData\n");
    fprintf(out,"2\n");
    fprintf(out,"\"%s \"\n",Name.c_str());
    fprintf(out,"\"DGInterpolationScheme\"\n");
    fprintf(out,"1\n");
    fprintf(out,"%g\n",t);
    fprintf(out,"4\n");
    fprintf(out,"%i\n",Step);
    fprintf(out,"%i\n",nbComps);                                                // nbEq //NOT YET IMPLEMENTED nb equations != 1
    fprintf(out,"%i\n",mesh->getNbElmt(eT));
    fprintf(out,"0\n");
    
    // --- Write solution ---
    
    int nbENP = (*fsIt)->getNbENP();
    
    for(int eID=0;eID<mesh->getNbElmt(eT);eID++){
      int elmtGmshID = mesh->getGmshTag(eT,eID);
      fwrite(&elmtGmshID,sizeof(int),1,out);
      fwrite(&nbENP,sizeof(int),1,out);
      for(int k=0;k<nbENP;k++){
        for(int dim=0;dim<nbVar;dim++){
          double u = (*sol)(eID,2);
          fwrite(&u,sizeof(double),1,out);
        }
        if(nbComps==3){double u = 0.; fwrite(&u,sizeof(double),1,out);}           // add 0 because vector is 3D!
      }
    }
    
    fprintf(out,"\"$EndElementNodeData\"\n");
    
  }  
  
  fclose(out);
}